Scaling exponents of Forced Polymer Translocation through a nano-pore 
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We investigate several scaling properties of a translocating homopolymer through a thin pore 
driven by an external field present inside the pore only using Langevin Dynamics (LD) simulations 
in three dimensions (3D). Motivated by several recent theoretical and numerical studies that are 
apparently at odds with each other, we determine the chain length (TV) dependence scaling exponents 
of the average translocation time (r), the average velocity of the center of mass (vcm), and the 
effective radius of gyration (R g ) during the translocation process defined as (r) ~ N a , (vcm) ~ N~ e , 
and R g ~ N" respectively, and the scaling exponent of the translocation coordinate (s-coordinate) 
as a function of the translocation time (s 2 (t)) ~ t 13 . We find a = 1.36 ± 0.01, f3 = 1.60 ± 0.01 for 

(s 2 (t) ~ t 13 ) and P = 1-44 ± 0.02 for (As 2 (r)> ~ t B , 8 = 0.81 ± 0.04, and v ~ v = 0.59 ± 0.01, 
where v is the equilibrium Flory exponent in 3D. Therefore, we find that (r) ~ TV 1 ' 36 is consistent 
with the estimate of (r) ~ (Rg) / (vcm) ■ However, as observed previously in MC calculations by 
Kantor and Kardar (Y. Kantor and M. Kardar, Phys. Rev. E, 69, 021806 (2004)) we also find 
the exponent a = 1.36 ± 0.01 < 1 + v. We also observe that a = 1.36 is in closer agreement with 
a — (1 + 2z/)/(l + v) as recently proposed by Vocks et al. (H. Vocks, D. Panja, G. T. Barkema, and 
R. C. Ball, J. Phys.: Condens. Matter 20, 095224 (2008)). We also discuss the dependence of the 
scaling exponents on the pore geometry for the range of N studied here. 
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PACS numbers: 87.15.A-, 87.15.H-, 36.20.-r 



INTRODUCTION 



Translocation of biopolymers accros a biomembrane, 
e.g., transport of RNA molecules out of a nucleus, in- 
vasion of viruses into cells, etc., are ubiquitous and 
important processes in biological systems [![. Recently 
voltage driven translocation of a single stranded DNA 
through a cv-hemolysin pore in biomembrane 0, and sub- 
sequently double stranded DNA through synthetic sili- 
con nanopores[3j have stimulated a lot of activities as 
the phenomenon is rich in fundamental science involved 
and its prospective technical applications for detecting 
DNA/RNA sequences. While it is the attributes of het- 
eropolymer translocation that are the key ingredients 
for prospective new sequencing methods, these exper- 
iments have generated stimulating theoretical and nu- 
merical studies directed towards fundamental physics of 
homopolymer translocation through a nanopore. An im- 
portant question that has been repeatedly raised is how 
does the average translocation time scale with the chain 
length and what is the equation of motion that describes 
the situation adequately. Approaches using Fokker- 
Planck equation with entropic barrier term incorporated 
in the free en ergy have generated useful insights to 
the problem More recently partial Fokker-Planck 
equation (PFPE) has been suggested is the natural lan- 



guage of the problem (il. [Tol| . Quite naturally, a number 
of simulational studies have been directed to test predic- 
tions of these theories 0- [13]. 

This paper is aimed at determining the relevant scal- 
ing exponents of forced translocation of a homopolymer 
through a nanopore by carrying out large scale Langevin 
dynamics (LD) simulations in three dimensions (3D) and 
comparing the findings with those predicted by theoret- 
ical arguments. We look at the arguments for the un- 
biased case first as it serves as the reference for extend- 
ing the theoretical arguments in the presence of an ex- 
ternal field. Naturally, the equilibrium radius of gyra- 
tion R g ~ N" of a chain of length N, where v is the 
Flory exponent, is used as the relevant length scale in all 
the theories. The first theoretical argument came from 
Chuang et al. [7( who predicted that for the unbiased 
translocation the mean translocation time should scale 
in the same manner as a freely diffusing chain so that 
(t) ~ R^/D ~ N 1+2v , assuming the diffusion coefficient 
D ~ 1 /N appropriate for the free-draining limit (no hy- 
drodynamic interaction) . In this theory it is argued that 
the Rouse relaxation serves as the lower bound and in 
presence of a nanopore a smaller amplitude accounts for 
the slowness of the process 1J| . This theory also pre- 
dicts that the scaling exponent of the reaction coordi- 
nate defined as (s 2 (r)) ~ t 13 is given by j3 = 2 /a. As 
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usualfj]- 11 1 we denote by s(r) the monomer that is in- 
side the pore at time r. Noticeably, the theory is essen- 
tially very simple and the exponents are functions of v 
only with a = 1 + 2v, = 2/(1 + 2v) so that a/3 = 2. 
In two dimensions (2D) v = 0.75 leads to a = 2.5 and 
/3 = 0.8 respectively. In three dimensions (3D) v = 0.588 
leads to a = 2.2 and (3 = 0.92 respectively. The the- 
ory put forward by Dubbledam et al. invokes an addi- 
tional surface exponent term t-| 25| so that for the dif- 
fusive case this theory predicts [9j, llOj ct = 2(1 + v) — 71 
and (3 — 2/ a. For unbiased translocation this theory also 
predicts the product a/3 = 2. Several recent numerical 
studies in 2D @, El, H, ljj and one in 3D[l8[ supports 
Chuang et al., while Dynamic Monte Carlo (DMC) re- 
sults by Dubbledam et al. report a — 2.5 and (3 = 0.8 in 
3D which contradicts Chuang et al. and supports their 
own prediction [9j. While all the simulation studies ver- 
ify a/3 = 2.0, recent theories by Panja et al. and Vocks 
et al. pointed out the role of decay time of monomer 
density near the pore and argues that the translocation 
time is anomalous up to the Rouse time £r ~ N 1+2v , and 
becomes diffusive afterwards [l2, 13 1. Therefore, for the 
unbiased translocation the collective numerical results do 
not support any of the proposed theories completely. 

Let us now look at the theoretical studies of driven 
translocation whose numerical verification including the 
underlying assumptions is the main focus of the paper. 
According to Kantor and Kardar (r) ~ (R g )/ (vcm) ~ 
N 1+u , assuming vcm ~ l/N. Kantor and Kardar[8j ar- 
gued that since the chain is only driven at one point inside 
the narrow pore, the accompanying change in its shape 
due to the bias is insignificant for the rest of the chain 
and therefore, the chain in this case is also described by 
the equilibrium Flory exponent v. To verify their scal- 
ing argument Kantor and Kardar carried out Lattice MC 
simulation of self-avoiding chains in 2D and noticed that 
the numerical exponent ~ 1.5 < 1 + v = 1.75. They 
argued that finite size effects are severe in this case and 
the relation (t) ~ N 1+l/ should be taken as an upper 
bound that will be seen only for the extremely large 
chains. Vocks et al. on the contrary, using arguments 
about memory effects in the monomer dynamics came 
up with an alternate estimate [3] (t) ~ N 1 +" . This 
seems to be consistent with most of the numerical data 
in 3D. However this estimate fail to capture the recent 
2D simulation results using Langevin dynamics and MC 
simulations [HI, [ijj where one sees a crossover of the 
a-exponent from 1.5 to 1.7 (as opposed to 1.428). Dub- 
bledam et al. have extended their PFPE based theory for 
the driven translocation[10] and came up with the follow- 
ing relations a — 2v + 1 — 71 and (3 = 4/ (2(1 + v) — 71). 
The prediction of Dubbledam et al. for the exponents 
are a = 1.55 and (3 = 1.56 in 2D and a = 1.5, and and 
(3 = 1.6 in 3D respectively. The DMC results of Dub- 
bledam et al. are consistent with this theory. However, 
more recent numerical results using LD and MD [22J, |23j 



produce similar results which are only in partial agree- 
ment with these theories. 

In this paper not only we calculate these scaling ex- 
ponents a and (3 for the driven chain but provide in- 
sights how the scaling aspects are affected by boundary 
and geometric factors by monitoring some of the rele- 
vant time dependent quantities during the translocation 
process. This allows us to check how well some of the as- 
sumptions are satisfied for the driven translocation and 
discuss possible scenarios for the disagreements between 
the theoretical predictions and numerical studies. Thus 
far these issues have not been adequately addressed in 
the literature. 



THE MODEL 



We have used the "Kremer-Grest" bead spring model 



to mimic a strand of DNA [24j . Excluded volume inter- 



action between monomers is modeled by a short range 
repulsive LJ potential 



U LJ (r) 
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Here, a is the diameter of a monomer, and e is the depth 
of the potential. The connectivity between neighboring 
monomers is modeled as a Finite Extension Nonlinear 
Elastic (FENE) spring with 

U FE NE(r) = -\kR 2 Q ln(l - r 2 /R 2 ) , 

where r is the distance between consecutive monomers, 
k is the spring constant and Ro is the maximum allowed 
separation between connected monomers. We use the 
Langevin dynamics with the equation of motion 

f^-Wi-rf - +Wi{t). 

Here T is the monomer friction coefficient and Wi(t), is 
a Gaussian white noise with zero mean at a temperature 
T, and satisfies the fluctuation-dissipation relation: 

< Wi(t) ■ Wj(t') >= 6fc B Tr% S(t - t') . 

The purely repulsive wall consists of one monolayer of 
LJ particles of diameter 1.5er on a triangular lattice at 
the xy plane at z = 0. The pore is created by removing 
the particle at the center. Inside the pore, the polymer 
beads experience a constant force F and a repulsive po- 
tential from the inside wall of the pore. The reduced 
units of length, time and temperature are chosen to be 
a, ery^", and e/ks respectively. For the spring potential 
we have chosen k = 30 and Rij — 1.5a, the friction co- 
efficient r = 1.0, and the temperature is kept at 1.5/fcs 
throughout the simulation. 



Q 



F=6.0 





H N=is6. F=6.0- 




■ N=128, F=6.0. 




■l N=256. F=4.0 




flL N=128, F=4.0" 


Lj 


/ i \tah..ii ' 




500 



1000 



1500 



FIG. 1: Histogram P(t) of flight times for chain lengths 
N = 64, 128, and 256 for bias F=6.0. The symbols repre- 
sent simulation data and the solid lines are fits with a form 
P(t) = At x exp(— /rr). The inset shows the corresponding 
scaled plots where the r-axis has been scaled by F(t) and 
the j/-axis has been scaled by the maximum value of the his- 
togram. 
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FIG. 2: Scaling of the mean translocation time (t) (logarith- 
mic scale) scaled by the applied bias F as a function of chain 
length TV (logarithmic scale). The open circles and squares 
refer to F = 4.0 and F = 6.0 respectively. The inset shows 
the corresponding scaling of vcm/F. 



We carried out simulations for chain lengths /V from 
8 — 256 for two choices of the biasing force F = 4 and 
6, respectively. Initially the first monomer of the chain 
is placed at the entry of the pore. Keeping the first 
monomer in its original position the rest of the chain 
is then equilibrated for times at least an amount propor- 
tional to the N 1+2v . The chain is then allowed to move 
through the pore driven by the field present inside the 
pore. When the last monomer exits the pore we stop 
the simulation and note the translocation time and then 
repeat the same for 5000 such trials. 

SIMULATION RESULTS AND THEIR 
INTERPRETATION 

Typical histograms for the passage time are shown in 
Fig. [1] for F — 6.0. When the time axis is scaled by 
the mean translocation time multiplied by the bias F 
and the peak of the distribution is normalized to unity, 
we observe (inset) a nice scaling of all the histograms 
on a single master curve. We also note that an excel- 
lent fit (solid lines) could be made with an expression 
P{t) — At x exp(— /at) for all the plots with the peak 
position being given by T max = A//i. We calculated 
the average translocation time from the weighted mean 
(t) = J Q maw tP(t)cLt, where t max for each distribution is 
chosen such that at t max the distribution P(t) is about 
0.01 % of its peak value. We have checked that (r) cal- 
culated from the area is marginally greater than T pea k 
obtained from P(t). 

The scaling exponent a of the mean translocation time 
(r) ~ N a is extracted by plotting the (r) as a func- 
tion of N shown in Fig. [2l Evidently, we find that 
(t) - l/F and (r) - TV 136 . The inset of Fig. EH shows 



that the velocity of the center of mass increases lin- 
early with the bias and scales as vcm ~ 1/iV ' 81 . We 
note that vcm does not scale as 1/N. It has been sug- 
gested that this exponent is not universal and depends 
on the width and the geometry of the pore[22j|. We will 
come back to this issue later. The scaling exponent (3 
of the s coordinate is shown in Fig. [3] For clarity, we 
have shown results for the two largest chain lengths N 
= 128 and 256. When we calculate the first and the 
second moments of s(r) we find that (s(r)) ~ t 8 and 
(s 2 (t)) ~ t 1 ' 6 for a wide range of the translocation time 
(the slope remains the same between the blue and the 
green vertical windows and between the green and the 
red vertical windows respectively in Fig. 3). The data 
as a function of the scaled translocation time Ft show 
excellent collapse. Since (s 2 (r)) ~ (( s ( t ))) 2 j one ex- 
pects to see (As 2 (t)} = (s 2 (t) — (,s(r)) 2 ) ~ r 1 - 6 during 
the same time window. However, (s 2 (r) — (s(r)) 2 ) re- 
veals additional features where the slope changes from 
(s 2 (t) — (s(r)) 2 ) ~ r 103 (between blue and green dashed 
vertical lines) to (s 2 (t) — (s(r)) 2 ) ~ t 144 (between green 
and red vertical lines). For the forced translocation 
(s(t)) ^ and it is likely that a tiny difference of time de- 
pendence of the first and 2nd moment during the translo- 
cation of the chain that is not visible in the plot of 1st or 
the 2nd moment of the s-coordinate becomes noticeable 
in its fluctuation. Therefore, if we use the fluctuations 
in s to define ((s(r) — (s(t))) 2 ) ~ r^, then from the late 
time slope (Fig. 3) then we get (3 = 1.44. 

We now compare these results with the theoreti- 
cal predictions and other existing numerical results. 
The translocation exponent in 3D, according to Kan- 
tor and KardarQ, is 1 + v = 1.588, and according to 
Dubbledam.0 is 1.50. First of all, as observed in 2D 
MC simulations by Kantor and KardarQ we also obtain 
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FIG. 3: variation of (s 2 (r)) (top, dotted) and (s(t)) (top, 
dashed-dot), and (s 2 (r) — (s) 2 ) (bottom) as a function of the 
scaled translocation time Ft. The black and blue colors cor- 
respond to chain length N=128 for F=4.0 and F=6.0 respec- 
tively. The red and magenta colors correspond to the chain 
length N=256 for F=4.0 and F=6.0 respectively. 



a smaller value of a = 1.36 ± 0.01 < 1 + v = 1.588 
in 3D. Kantor and Kardar argued that a lower value is 
obtained due to finite size effects and expect that for 
very large chains one would find 1 + v ~ 1.59. This 



bound has recently been criticized by Vocks et al. [13J] who 
using arguments about memory effects in the monomer 
dynamics came up with an alternate exponent estimate 
given by = 1.37. Evidently, our result is in agree- 

ment with this prediction. As for the exponent (3 we find 



(s 2 (r)> ~ t 5 and (As 2 (r)) = (s 2 (r) - (s) 2 ) 



(if 



we use the later window). Therefore, with Dubbledam 
et al. we do not agree with the calculated value of a, 
but Dubbledam et al. also used (s 2 (t) ~ t' 3 , to de- 
fine the exponent f3 and the numerical value /3 = 1.6 
is exactly the same as found here. It is noteworthy 
that the fluctuation (As 2 (r)) is time dependent and the 
slope of (As 2 (t)} ~ r 103 at early time crosses over to 



<A s 2 (r)} 



1.44 



at a later time, while the slope for 



(s (t)} ~ t is constant for a wider range. If we use 
(3 = 1.44, obtained from the definition of fluctuation of 



the s coordinate, then we find the relation a/3 — 2.0 is 
satisfied for the forced translocation as well. This trend 
is qualitatively the same for the simulation using a square 



pore[23j, where we find that (r) ~ A 1 



,-1.52 



and (As 2 (r)) ~ r 145 (so that a/3 ~ 2.0, same as re- 
ported here if we extract j3 from the slope of the plot 
(As 2 (r)) ~ r). Our results may be relevant in the con- 
text of a recent recent article by Chatelain, Kantor, and 
Kardar [27] who showed that the variance of the proba- 
bility distribution P(s,t) grows subdiffusievly. 

We now look more closely at the factors responsi- 
ble for the translocation process. The expression r ~ 
(R g ) I (vcm) ~ N 1+u has two components: the depen- 
dence of Vcm on N and R g on N respectively. We now 
look at these two components separately. During the 
driven translocation the chain does not find enough time 
to relax. Therefore, it is important to know how does 
the shape of the chain vary as a function of time and 
how different it is compared to its equilibrium configura- 
tion. During the forced translocation at any instant of 
time only one segment of the entire chain feels the bias. 
Kantor and Kardar Q argued that the shape of the chain 
is hardly affected by it so that it will still be described 
by the equilibrium Flory exponent v. This argument will 
not be strictly valid for the model used here as the beads 
are connected by elastic bonds and it is expected that 
quite a few neighbors on either side of the driven bead 
inside the pore will be indirectly affected by it. 

In order to verify this issue first, we have calculated 
the equilibrium (R g ) of the chain clamped at one end 
at the pore in presence of the same LJ wall. We 
find v ~ 0.6 ± 0.01 (Fig. SJ). We have also calcu- 
lated the relaxation time r r of the end-to-end vector 
((Ri N (t + r) - (R)) • (R 1N (t + r) - (R))) ~ exp(-t/r r ) 
and checked that we get the same v from the relaxation 
measurements. This is consistent with the theoretical 
prediction of Eiscnrieglcr, Kremer, and Binder that in 
presence of the wall the exponent v remains the same as 
that of its bulk counterpart 26] . To get an idea how fast 
is the translocation process, compared to the correspond- 
ing relaxation time, for the chain lengths N = 64, 128, 
and 256, we find r r - 1000,4500, and 20200 respectively 
and the corresponding average translocation times (r) 
are 215, 530, and 1330, respectively. Even in the linear 
response regime where (r) ~ 1/-F 1 , we observe r r >> (r). 
The insets of Fig. |4] shows the time dependence of R g (t) 
(we use a different notation R g for the driven chain) . We 
notice that during the translocation process the chain is 
significantly elongated around t ~ 0.5 (r) and acquires 
relatively compact structure immediately upon exiting 
the pore. The dashed lines (black and green) show the 
corresponding average values R g from which we extract 
the exponent v ~ v (Fig. [4|). Contrary to what is as- 
sumed by Kantor and Kardar, we notice significant dis- 
tortion of the chain. Surprisingly, we find that the R g 
scales almost the same way as (R g ) ~ A 6 (slopes are 
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FIG. 4: Equilibrium R g and effective R g during the translo- 
cation process. The absolute value of the effective R g is larger 
than the equilibrium value as it is pulled, but both of them 
scale with the equilibrium Flory exponent. The insets shows 
the average time dependence of the woiwr(t) and R g (t) during 
the translocation. The straight lines in the inset represent the 
average value. 



the same in Fig. 4). even when (r) « r r . Therefore, 
numerically we find that the chain is still described by 
the equilibrium (R g ). 

Likewise, as expected in LD simulation, we notice 
that vcNi(t) saturates quite quickly and this value is 
almost the same during the translocation process and 
~ (v CM ) ~ l/N 5 . where 5 = 0.81 ± 0.04. Since v ~ v, 
our studies indicate that it is the (vcm) which does not 
exhibit inverse linear dependence on chain length N is the 
responsible factor for the deviation from (r) ~ N 1+v . 

It is worth mentioning that we have carried out exactly 
the same LD simulations with wall particles on a square 
lattice^ . We find that (r ~ N 1A \ (s 2 (t)} ~ t 1 - 52 , 
and (As 2 (t)) ~ t 1,45 (so that a(3 ~ 2.0, same as re- 
ported here if we extract (3 from the slope of the plot 
(As 2 (t)} ~ r). These results for the square pore are also 
consistent with recently reported LD and MD simulation 
results in 3D using GROMACSjH]. Recently Gauthier 
et al. [2l| carried out similar studies of polymer transloca- 
tion through a narrow pore (including hydrodynamic in- 
teractions) and found a systematic variation of the mea- 
sured scaling exponents as a function of the pore width. 
However, their studies are limited to relatively narrow 
range of N up to 31 only. In our studies the exponents 
for a relatively wide range of N seem to depend on the 
pore geometry. Whether this implies true nonuniversal- 
ity or not remains an open issue. 



CONCLUSION 

To summarize, we have used Langevin dynamics in 3D 
to study the scaling properties of a driven translocating 
chain through a nanopore. We notice that the chain un- 



dergoes a significant shape change during the fast translo- 
cation process, contrary to what assumed by Kantor and 
Kardar is formulating the theory of forced translocation. 
However, despite significant distortion, we find the chain 
is still described by the equilibrium Flory exponent. We 
find that the (vcm) does not scale as its bulk counter- 
part and depends on pore width and geometry. It is 
likely that density variation on either side of the pore 
during the translocation process affects the overall mo- 
tion of the chain. We find that the a = 1.36 < 1 + It is 
worth mentioning that the collective numerical work by 
various groups failed to validate the Kantor and Kardar 
result a = 1 + v for the forced translocation, including 
the results listed here. Likewise, although the value of 
a = 1.36 that we obtain is in excellent agreement with 
the analytical estimate of Vocks et al. a — 37^7 = 1-37 
in 3D, the results from 2D simulations do not agree with 
the estimate of Vocks et al.. Finally, we notice a differ- 
ence in the s-exponent (3 when calculated from its second 
moment ((3 = 1.6) and its fluctuations ((3 — 1.44). The 
later (J3 = 1.44) agrees with a/3 ~ 2.0 while (3 = 1.6 
overestimates it (a/3 ~ 2.2 > 2.0). The fluctuations 
(As 2 (t)) seem to reveal more structures not adequately 
studied so far. When we compare these results with the 
existing theories and other numerical results we notice 
that these results only partially support one theory or 
the other. Certainly more numerical and analytic work 
are needed for a more comprehensive understanding of 
forced translocation through nanopore. 
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